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We investigate theoretically the nonlinear state of ideal straight rolls in the Rayleigh- 
Benard system of a fluid layer heated from below with a porous medium using a Galerkin 
method. Applying the Oberbeck-Boussinesq approximation, binary mixtures with pos- 
itive separation ratio are studied and compared to one-component fluids. Our results 
for the structural properties of roll convection resemble qualitatively the situation in the 
Rayleigh-Benard system without porous medium except for the fact that the streamlines 
of binary mixtures are deformed in the so-called Soret regime. The deformation of the 
streamlines is explained by means of the Darcy equation which is used to describe the 
transport of momentum. In addition to the properties of the rolls, their stability against 
arbitrary infinitesimal perturbations is investigated. We compute stability balloons for 
the pure fluid case as well as for a wide parameter range of Lewis numbers and separation 
ratios which are typical for binary gas and fluid mixtures. The stability regions of rolls 
are found to be restricted by a crossroll, a zigzag and a new type of oscillatory instability 
mechanism, which can be related to the crossroll mechanism. 



1. Introduction 

The Rayleigh-Benard system is one of the classical setups and possibly the most pop- 
ular one to study pattern formation and hydrodynamic instabilities under well-controlled 
conditions. The system consists of a fluid layer bounded by two plates and is heated from 
below. If the temperature difference between the plates exceeds a critical value, the fluid 
leaves its static state and thermal convection sets in. Under appropriate conditions, the 
fluid flow forms a regular pattern, for example a state of ideal straight rolls (ISR), also 
called Benard rolls or steady overtuning convection (SOC). 

Over the last decades, the roll convection and in particular the instability mechanisms, 
which tend to limit the s tability region of the ideal roll s tate, have been investigated i n 
detail experimentally (see lBusse fc Whitehead! (Il97ll) andlCroquette fc Williams! (Il989h) 
as wel l as theoreticall y (sec Buss e fc Clever (ll979l):ISchluter et all (ll965l) : IClever fc Bussd 
(|l990l ): iBolton et all (|l986h : iBolton fc Bussd fjl985h ). One main result of this research 
is that the region of rolls, which are stable against infinitesimal perturbations, is re- 
stricted by certain instability mechanisms, for instance the zigzag, the skewed varicose, 
the Eckhaus and the crossroll mechanism. This remains also valid if the fluid between 
both the plates sat urates a por ous medium, a case apply ing to many natural and indus- 
trial proc esses (see Vafai ( 2005) or Nield fc Beiarj ( 20061 ) for an overview). For a porous 
medium, De La Torre Juarez fc Bussd ( 19951 ) investigated theoretically the stability of 
rolls against two-dimensional perturbations and found the Eckhaus and an oscillatory 
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instability mechanism. Straus! ( 19741 ) did the same but allowed for three-dimensional 
perturbations and found a zigzag as well as a crossroll mechanism. These investigations 
we re complemented by t he ex perimental observa tion of roll convection in porous media 
bv lShattuck et all (|l997i ) and lHowle et all (|l997l ). 

The aforementioned research was carried out for a one-component, i.e., a pure fluid. 
However, numerous applications are known where the fluid has to be treated as a 
binary mixture consisting of two different fluids (see for example Platten 20061) . Es- 
pecially, if the binary mixture shows a non- vanishing Soret effect, i.e., if concentra- 
tion currents are driven by temperature gradients, the dynamics of the system are ex- 
pected to change due to a coupling of the temperature field into the concentration field. 
For the classical Rayleigh-Benard system without porous medium, there has been a 
large body of work d ealing with binary m i xtures under th e influence of the Soret effec t 



Eaton et all (Il99lh: iBarten stall dl995h 
Dom inguez-Lerma et al. 
198 



Huke. fc Luckel (l2002h: iLe Gal et all (Il985h; 
1995D:lBarten et al. dl989l):ISchopf fc Zimmermannl (ll993h:lKnobloch fc Moore 
Cross fe Kiml(ll988l ) ;lKolodner et all ll986);lFutterer fe Luckel (|2002fklT ouiri et all 
Ahlers fc Rehberd (|l986l ): IWalden et all (|l985h ). In particu lar, the roll convec- 



tion and the corresponding stability mechanisms are well known (see Huke et al. ( 2000l )). 
For the system with porous medium, however, the standard of knowledge is much less 
developed, although there has been some work conce r ning t he stability of t he g round 
state and monocellu lar flow bv lCharrier-Moitabi et al. ( 2007 ). [Elhaiiar et al. (|2008l ) and 



Sovran et all (|200ll ). 

With this paper, we aim at extending the knowledge about flow patterns of binary 
mixtures in porous media by analysing the structural properties as well as the stability 
of the Benard rolls. In doing this, we restrict ourselves to a positive Soret coupling in 
which the Soret effect destabilises the ground state. The paper is organised as follows: 
In Sec. II we will briefly explain the basic equations of the system, the Galerkin method 
that we have used for numerical computations, the ground state and its stability. In Sec. 
Ill we discuss the properties of roll convection by investigating the structure of the fields, 
the advective heat transport and the mixing. Sec. IV contains the main results of our 
paper. Here, we represent the stability boundaries of the Benard rolls for a pure fluid as 
well as for binary mixtures and compare the former ones to the known results. In Sec. V 
we conclude with a summary of our results. 



2. Foundations 

2.1. System and basic equations 

We consider a horizontal layer of a porous medium filled with a pure fluid or a binary 
mixture in a homogeneous gravitational field, g = — ge z . The layer has thickness d and 
a vertical temperature gradient is imposed by fixing the temperature 

T = T ±— at z = T^ (2- 1 ) 

which can be realized in experiments by using highly conducting plates. T is the mean 
temperature of the layer and we assume AT > 0, i.e., the lower plate has a higher tem- 
perature than the upper one. Moreover, we consider the plates to be infinitely extended, 
rigid and impermeable. The porous medium is treated as isotropic and homogeneous. 
Furthermore, we assume local thermal equilibrium between the fluid and the porous 
medium such that there is no heat transfer between both phases. 

Convection is characterised by the fields of temperature T, the Darcy velocity v (also 
called seepage velocity), mass concentration C of the lighter component, total mass den- 
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sity p = pi + p2 and pressure P. Applying the Oberbeck-Boussinesq approximation, we 
assume the dynamic viscosity, the thermal expansion coefficients and the heat capaci- 
ties to be constant and equal to their values at the spatial averages To, Co, Pq of the 
thermodynamic variables. Moreover we neglect effects like radiation or the generating of 
heat due to friction between the fluid and the porous medium. For binary mixtures, the 
Soret effect can play an important role. We incorporate this effect in our model, which 
describes the generation of concentration currents due to temperature variat ions. Then, 
the balance equations for our system read according to lNield fe Beianl (|2006f ): 



V-v = (2.2a) 
c aPo d t v = -VP -^v + p [l + p T T + p c C]ge z (2.26) 
C tot d t T + C f {v-V)T = \ tot X7 2 T (2.2c) 

<f>d t C + (v • V)C = AotV 2 C + Aot^V 2 T, (2.2d) 

Jo 

Here, po is the mean density of the fluid, rj the dynamic viscosity, f3r (Pc) the thermal 
(solutal) expansion coefficient of the fluid and K the permeability of the porous medium. 
C denotes the heat capacity per unit volume according to its subscript /, s or tot of 
the fluid, the solid matrix or of the total medium (fluid and solid matrix). The heat 
capacities per unit volume can be connected via the porosity <f> of the porous medium, by 
the relation C to t = 4>Cf + (1 — <j>)C s . The same holds true for the thermal conductivities 
Xf , A s and Atot • -Dtot is the concentration diffusivity of the total medium and is equal to 
cj)D whereas D is the concentration diffusivity of the fluid. Jct is the thermodiffusion ratio 
which characterises the strength of the Soret effect. The correction factor c a emerges in 
front of dt\i to brin g the D arcy equation (|2.2b ) in agreement with experimental results 
(see iNield fe Beianl (|2006h ). Since the porous medium is assumed to be isotropic and 



homogeneous, c a is a scalar. 

In the basic equations (|2.2p we scale lengths by d, time by the vertical diffusion time 
^gf , temperature by Cf ^ gpo , concentration by Cf ^ gpo and pressure by . 
Moreover, we introduce the reduced deviations of the velocity u = (u, v, w), temperature 
9, concentration c and pressure field p from the conductive state thus obtaining the 
following set of balance equations: 

V-u = (2.3a) 

la d t u = -Vp - u + (9 + c)e z (2.36) 

d t 9 + (u • V)0 = Rw + V 2 9 (2.3c) 

4>*d t c + (u • V)c = Rtjiw + L (V 2 c - VV 2 6») . (2.3d) 

The normalised porosity (j>* = (f)Cf/C to t as well as the correction factor j a = JicTlpo ° a 
are scalars for an isotropic and homogeneous porous medium, (jf can have values in t he 
interval [0, 1] whereas 7 a is very small for usual porous materials ( Nield fc Beian ( 20061 )). 



Because of the smallness of 7„ , the time-derivative term in the momentum eq uation is of- 



ten neglected (see for example ISovran et all (120011) and Elhaiia r et all (120081) ). However 



we ret ain this term since it is - according to IVadasz fc Olekl (| 1999T) or IVadasz fc Olek 



( 20001 ) - essential for certain stability analyses. Note that (jf also appears only in front of 



a time derivative. Therefore, the correction factor and the normalised porosity can play 
a role for time-dependent phenomena only. 
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Via the concentration diffusivity Z?tot, the Lewis number L is defined by 

L= D^C L ^ (24) 

and therefore compares the time scales of concentration and heat diffusion. Note that 
the Lewis number as defined here is the inverse of the Lewis number Le as defined by 
Nield fc Beianl ()2006h . The Rayleigh-Darcy number is given by 



R = Po^rKC f d AT 
and measures the thermal driving. The separation ratio 

i> = -4^ k T (2-6) 

is proportional to the thermodiffusion ratio kx thus incorporating the Soret effect which 
emerges via —Li[;V 2 9 and the term Ripw in the concentration balance equation (|2.3t L. 
In this paper, we cover mainly separation ratios from to 0.5, which are typical for 
alcohol-water mixtures at high alcohol conc entration (30-80%) a nd many gas mixtures, 
like Ne-C0 2 , He-Xe, H 2 -Xe and Ne-Ar (see iLiu fc Ahlersl (|l997l )). That means in par- 



ticular that we restrict ourselves to a positive Soret coupling {ip > 0), where the lighter 
component of the mixture is driven into the direction of higher temperature. 

In the momentum equation (|2.3b ). which is deduced from Darcy's law, advective 
transport and diffusion of momentum are neglected. Instead, we take into account the 
relaxation term — u describing the friction between the fluid and the porous matrix, 
which dominates the aforementioned mechanisms at low Reynolds numbers (Re = 0(1)). 
Thereby, to define the Reynolds number of a flow through a porous medium, the average 
radius of the pores is chosen as the characteristic length scale. 

Since the fluid is assumed to be incompressible, we can write the velocity field u as 

u = V x V x $e 2 + V x *e 2 , (2.7) 

which automatically fulfils the mass balance (|2.3a ) . To derive equations for the potentials 
$ and 'J, one can apply the curl to the momentum balance equation p. 3b ) once or twice, 
respectively and take into account only the third component of the obtained equations. 
This also eliminates the pressure term. To derive boundary conditions for we assume 
that there is no external pressure gradient which causes a meanflow U(z, t) — (u(x 1 z, t)) x , 
with being the spatial average in x-direction. Averaging (|2.2b ) over x shows that 
any meanflow will decay on the time scale 7 a . The same can be demonstrated for "J. Not 
being interested in fast transients, we assume that U = ^ = 0. 

At the impermeable plates, the vertical concentration current vanishes, i.e., d z (c—ip9) = 
for ±i. To avoid this coupling of the fields in the boundary conditions, we introduce in 
a next step the field 

( = c~il>e. (2.8) 

The corresponding balance equation for £ can directly be derived by combining the 
equations (12.3b ) and (|2.3H ). Altogether, the new balance equations read 

~f a d t V 2 A xy <5> = -V 2 A xy $ ~ A xy {9 + c) (2.9a) 
d t + (u • V)9 = Rw + V 2 9 (2.96) 
4>*dtC + (u • V) C + (1 - <P*) (u ■ V) 9 = (cj)* - 1) RtpA xy ^ + LV 2 C - 0>V 2 6» (2.9c) 
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with As y := &l + dy and the boundary conditions 

= $ = (9 = <9 Z C at z = ±-. (2.10) 

2.2. Numerical method 

To obtain roll solutions, we used a Galerkin method with the following ansatz for the 
fields X = $,6, C 

X(x,z) = X mn cos(mkx) f n (z) . (2-11) 



This ansatz is almost the same as in iHuke et all (j2000h where it is described in more 
detail. 

The difference is that ip = and that <E> is expanded in the vertical direction by 

* f ( \ Jv 7 ^ cos(nTrz) if n = l,3,5... 

I v 2 smynitz) if n = 2,4,6... 

to satisfy the boundary condition (|2 . 10|) . Rolls are even in x with an approp riate choice of 
the origin and fulfil the so-called mirror-glide symmetry (see Veronisl ( 19661 ) ). This allows 



to simplify the ansatz (|2.1ip by dropping half the modes. In the following computations, 
all 6- and £-modes with n + m > N and all <I>-modes with n + m > N/2 are neglected, 
where we choose iV = 24 if not otherwise specified. 

The stability of rolls against arbitrary infinitesimal perturbations is tested by expanding 
the perturbations as follows 

6X(x,y,z,t) =J2 SX ™ne st+i[ ( d - mk)x+by] f n (z), (2.13) 

m,n 

where the ansatz is truncated in a way consistent with the truncation of (|2.11[) . Except 
for the differences pointed out already, this ansatz and the corresponding method can be 
found in Huke et al. (2000). The symmetries of the underlying roll pattern allow for 6 = 



to divide the perturbations into perturbations that are even or odd in the x-direction 
and into G- as well as G-perturbations, that reduce to symmetric and antisymmetric 
perturbations under the mirror glide operation in the special case d=b=0. Moreover, 
because of these symmetries, only perturbations with b > and d £ [0, k/2) have to be 
tested. 

2.3. Ground state and linear Stability Analysis 

The linear stability problem f or binary mix t ures o f the system under consideratio n has 
already been i nvesti gated by ISovran et al. ( 2001 ). Charrier-Moitabi et al. ( 2007 ) and 



Elhaiiar et al. (2008). Having used an analytical method and a Galerkin expansion, we 



obtain results in good agreement with the former. In the following, we are going to 
review the most important facts concerning the conductive state as a preparation and 
introduction for the subsequent nonlinear analysis. 

In the ground state, the fluid rests and heat is transported only by diffusion. The 
temperature difference AT between the plates imposes a linear temperature profile as 
follows: 

AT 

Tco„d(z) = To - — z. (2.14) 
Due to the Soret effect, this temperature gradient generates a linear concentration profile 
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of the form: 



C oooA (z) = C + ^£fz. (2.15) 



The pressure distribution in the conductive state is given by 
Pcond(z) =Po + P9 



(2.16) 



which can be calculated from the concentration and temperature held via the momentum 
equation (|2.3b ). If there is a positive Soret coupling as investigated here, the lighter 
component of the binary fluid is driven to the lower plate. This increases the density 
difference between the plates and therefore destabilises the ground state. 

The destabilising effect is illustrated in Figure Q] where the critical Rayleigh-Darcy 
number R c as well as the corresponding critical wavenumber k c is plotted against the 
separation ratio -0 for several L. For the pure fluid (ip = 0), the ground state loses its 
stability above = Air 2 w 39.48 against stationary perturbations of a lateral wavenum- 
ber = it. By contrast, the ground state of binary mixtures becomes already unstable 
for a thermal driving weaker than R® and the critical Rayleigh-Darcy number decreases 
with stronger Soret effect, i.e., growing ip. In particular, if the concentration gradients 
are slowly diffused away, i.e., for small Lewis numbers, the destabilisation is especially 
strong. The corresponding critical perturbations remain stationary but their wavelength 
k c goes to zero for large if). For example at L = 0.5, one finds k c — for ab out i/> ^ 1.759 
which is close to the theoretical value ip ^ atr = 1.7586 obtained by ISovran et al 



(2001). In experiments, the critical wavelength will then be as large as the finite size of 



the convection cell allows. 



3. Properties of the roll convection 

3.1. Structure of the fields 

To understand roll convection of a binary mixture from a qualitative point of view, 
one can study the temperature and concentration distribution as well as the flow field. In 
what follows, we give such a qualitative description by Figure[2] in which the streamlines, 
the concentration and temperature field along with the lateral profiles of the fields at 
midheight z = are shown for several values of thermal driving. Since the behaviour 
remains in principle the same over a wide range of wavenumbers, we fix the wavenumber 
to k — k® — it. The chosen parameters, L = 0.01, -0 = 0.3, can be realized easily 
in experiments with alcohol-water mixtures. In this case, the critical Rayleigh-Darcy 
number R c is reduced to about 0.39. The left column in Figure [2] displays the fields 
over one periodicity interval in x-direction for R = 10. This value lies in the so-called 
Soret regime where convection is dominated by the Soret effect. The temperature field 
deviates only marginally from its linear profile in the ground state, i.e., advective heat 
transport is weak. In contrast, the concentration field is already strongly modulated and 
forms plume- like structures since the slow concentration diffusion (L <C 1) allows for a 
perturbation of the linear profile already for weak advection. The anharmonicity of the 
concentration field becomes obvious by looking at the lower part of Figure [H where the 
lateral profile of the concentration field is plotted versus x at midheight, z = 0. The 
velocity field, represented by its z-component at midheight, is also anharmonic and the 
streamlines that illustrate the roll-like flow are deformed. 

The middle column of Figure[2]refers to R = 50. This value lies in the so-called Rayleigh 
regime where convection would set in also without an operating Soret effect. Here, a 
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stronger modulation of the temperature field can be seen due to an increased advective 
heat transport. Nevertheless, the temperature field is almost harmonic as seen from the 
horizontal variation of T at midheight. The concentration field shows the characteristic 
boundary layer behaviour: The binary fluid is well mixed in the bulk whereas pronounced 
concentration gradients exist at the plates and the roll boundaries. The better mixing 
is caused by the larger velocity of the fluid. The corresponding streamlines reflect the 
roll-like flow and are almost harmonic indicating that the anharmonic behaviour of the 
velocity for R = 1 is induced by the Soret effect. This conclusion has been strengthened 
when we simulated the streamlines of a pure fluid and did not find any deformation 
provided the thermal driving was not too strong. Note furthermore that the deformation 
of the streamlines in the Soret regime seems to be generic for the system with porous 
medium. This can be explained by the replacement of the momentum diffusion term V 2 u 
(contained in the Navier-Stokes equations) by the relaxation term — u in our momentum 
balance equation (|2.3t >). In the clear fluid, i.e., in the system without porous medium, 
momentum diffusion tends to smooth the spatial anharmonicity of the velocity field 
caused by the very anharmonic concentration field. Since the relaxation term — u in 
(|2.3t >) does not contain spatial derivatives, it does not provide spatial smoothing and 
the anharmonicity of the concentration field is imposed via the buoyancy term onto the 
velocity field. Taking into account that the buoyancy at R = 50 is caused rather by the 
almost harmonic temperature field than by the concentration field, it becomes clear why 
the streamlines are barely deformed. 

When nonlinear effects are amplified further by increasing the heating rate, plume-like 
structures appear also in the temperature field and the concentration boundary layers 
become thinner. We can see this in the right column of Figure [21 where R has reached a 
value of 100. The streamlines are now deformed again, this time because the temperature 
field has become anharmonic by the intensive thermal driving. 

3.2. Nusselt and mixing numbers 
To describe the advective heat transport due to roll convection quantitatively, we make 
use of the Nusselt number Nu. The Nusselt number is given by 

Nu . = foot,,) (3 1} 

Vcond,2/ 

where jtot,z denotes the total vertical heat current density and j C ond,z the vertical heat 
current density in the conductive state. (•) denotes the lateral average. The Nusselt 
number at the plates can be computed from the modes obtained by our Galerkin method 
with truncation index N as follows: 

F)tt N 

Nu=l-^^(-l)"n# 02 „. (3.2) 

n— 1 

Note that the Nusselt number does not actually depend on the z-position in a stationary 
state of convection as the rolls are. 

In Figure [3] we compare Nusselt numb ers of the pure fluid, which w e have found to be 



in good agreement with the results from lDe La Torre Juarez fc Bussd (|I 9951 ) . to Nusselt 



numbers of binary mixtures. In the ground state, Nu is equal to I, since heat is trans- 
ported only by diffusion. After the roll convection has started, i.e., R has exceeded R c , 
the Nusselt number and simultaneously the advective heat transport increases mono- 
tonically with R. For a given Rayleigh-Darcy number, Nu is always larger for binary 
mixtures with ip > than for the pure fluid, as the Soret effect causes a concentration 
gradient giving rise to a more pronounced buoyancy and thus a stronger convection. For 
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the same reason, the Nusselt number grows when the Soret effect becomes As the advec- 
tive mixing is strong in the Rayleigh region, the Nusselt number of a binary mixture with 
small L approaches the one of the pure fluid there. However, when L becomes larger, i.e., 
when concentration diffusion is fast, the advective mixing does not succeed to mix away 
the concentration gradient and the differences in Nu to the pure fluid case remain more 
significant. 

We have also found these tendencies for wavenumbers that differ from the wavenumber 
k = 7r chosen in Figure albeit convection becomes weaker and Nu decreases if the 
wavenumber gets too large or too small. 

The mixing of a binary fluid can be described by the so-called mixing number M, 
which is defined by the normalised variance of the concentration field: 



„ ffif (33) 

V (^cond) - (Ccond) 2 

C^ond denotes the concentration field in the conductive state and (•) the spatial average. 
From M = 1, which is the value of the mixing number in the ground state by definition, 
it decreases when roll convection sets in and the components of the binary fluid are 
mixed. According to Figure [H the mixing number decreases only slightly in the Rayleigh 
regime when the thermal driving gets stronger. In other words, one cannot mix the 
fluid components perfectly by increasing R. The reason is that the fluid is already well 
mixed in the Rayleigh regime except for the fluid layers near to the plates and at the roll 
boundaries, whose thickness depends barely on R. Instead, the boundary layer behaviour 
in the Rayleigh region depends mainly on the Lewis number of the mixture such that the 
mixing number is reduced when L is lowered. The latter facts are also valid in the clear 
flu id case and can be explained using concentration boundary layer theory in analogy 
to Hollinger (|l996h . A larger separation ratio improves the mixing mainly in the Soret 



regime whereas the influence of ip in the Rayleigh regime is rather weak. 



4. Stability of the rolls 

4.1. Stability boundaries for the pure fluid 

Using the Galerkin method from section [2721 with N = 24, we have tested the roll struc- 
ture against infinitesimal perturbations of arbitra ry wavenumber. For the pure fluid , 
simi lar stab i lity a nalyses have been carried out by Dc La Torre Juarez fc Bussej (|l995h 



and Straus! (1974). In agreement with their results, we find that only two instability 



mechanisms, the zigzag and the crossroll mechanism, limit the region of stable rolls. The 
corresponding stability boundaries in the (R, fc)-parameter space are shown in Figure [5j 
Below, we review briefly the important instability mechanisms and compare our results 
to those of Straus and Juarez. If not otherwise specified, the instabilities are related to a 
real eigenvalue, i.e., the corresponding stability boundary does not depend on <jf or j a . 

4.1.1. Zigzag instability 

The zigzag (ZZ) boundary, denoted by the dash-dotted line in Figure restricts the 
region of stable rolls on the small-/c side; if the rolls are zigzag-unstable, a new set of rolls 
with a larger wavenumber begins to grow. The ZZ perturbations belong to the subclass 
of G-perturbations that are odd in x-direction. They have the same periodicity in re- 
direction as the existing pattern. Consequently, we find the ZZ instability for d = 0. It is 
sufficient to check the stability only for a single point on the b-axis near b = in order 
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to det ermine whether a roll state is ZZ-stable or not (for further details see lBolton et 
(|l98rf l). 



Compared to the ZZ boundary calculated by Straus, ours is much more restrictive. At 
R = 55 for example, rolls should be unstable for k < 3.046 according to our computations 
whereas Straus states the same for k < 2.5 using an analytical criterion. 

4.1.2. Crossroll instability 

According to Figure [5l the region of stable rolls is mainly restricted by the crossroll 
(CR) boundary given by the solid line. Thus, rolls with either too large or too small 
wavenumbers are destabilised by the CR instability, which causes the growth of rolls per- 
pendicular to the existing pattern. Above R ss 342 =: fi2, the CR mechanism destabilises 
the rolls independent of their wavenumber. The CR perturbations are even in x-direction 
and belong to the subclass of G-perturbations. They have the same periodicity in the 
IE-direction as the original rolls so that we find the CR instability for d = 0. By contrast, 
the parameter &, which represents the wavenumber of the CR perturbation perpendic- 
ular to the roll pattern, cannot be fixed to test wether rolls are CR-stable. Instead, to 
decide this, one has to find the value b max where the most critical eigenvalue reaches its 
maximum by applying an interpolation procedure. In Figure the value of b max on the 
CR boundary is plotted as function of k, whereby the upper section of b max (k) belongs 
to the part of the CR boundary at higher R. Near the onset, b max is close to fc°. For 
larger R, b max is generally larger since the whole stability region shifts to larger k. 

The above mentioned value R° h « 342 differs from the one obtained by Straus who 
states w 380. Furthermore, our CR boundary is again more restrictive on the small-fc 
side: According to our calculations, rolls with a wavenumber k < 2.28 are CR-unstable 
whereas Straus finds CR-stable rolls with few 1.8. However, these disagreements decrease 
when we lower our truncation parameter iV and thus we conclude that our results are 
more precise than those obtained by Straus. 

4.1.3. Eckhaus instability 

The Eckhaus (EC) boundary, denoted by the dotted line in Figure [5j lies below the 
CR boundary for all wavenumbers, i.e., the EC perturbations can only grow where the 
rolls are already CR-unstable. The EC instability tends to establish rolls with a better 
wavenumber in the direction of the wavevector of the original roll state. The correspond- 
ing perturbations are even in x-direction and fall into the subclass of G-perturbations. 
As purely two-dimensional perturbations, they are found at b = 0. Conveniently, the 
question of EC stability can be answered by i nves tigating a single point on the d-axis 
near d — (see De La Torre Juarez fc Bussd (1995)). Our result for the EC boundary 



agrees well with the one obtained by Juarez. 

Since we have found no other instability mechanisms inside the ZZ- and CR-stable 
region, we conclude that within this stability region rolls are stable against arbitrary in- 
finitesimal perturbations. In particular, rolls are stable at the critical point (i?°, fc°). Note 
that the stability region discussed above is similar to the stability ballon of the clear fluid 
at larg e Prandtl numbers for No-slip boundary conditions comp uted bvlBusse fc Cleverl 
(1979) as well as for Free-slip boundary conditions computed by lBolton fc Bussd (|l985f h 



In each of the three cases, the region of stable rolls is restricted by the ZZ boundary 
on the small-fc side and otherwise by the CR boundary. The similarity of the stability 
balloons can be understood using a heuristic argument. The influence of the advective 
term and the time derivative term in the Navier-Stokes equation becomes small for large 
Prandtl numbers so that it resembles the Darcy equation for small j a except for the fact 
that the relaxation term in (|2.3b ) is replaced by a diffusion term. 
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4.2. Stability boundaries for binary mixtures 

Figure [6] shows the EC, ZZ and CR boundary for binary mixtures below R — 80. We cover 
a range of separation ratios ranging from 0.01 to 0.4 and Lewis numbers from 0.1 to 1. 
Again, the stability boundaries are found to be in qu alitative agr e ement with those of the 
clear fluid for large Prandtl numbers, computed by Huke et all (2000): In the Rayleigh 



regime, the ZZ boundary Rzz(k) lies close to the corresponding line of the pure fluid and 
shifts slightly to smaller k if if> is increased. At the transition between the Rayleigh and 
the Soret region, Rzz(k) bends towards small k (see for instance L = 0.5 and if) = 0.05). 
By decreasing R further on, the ZZ boundary terminates with a finite slope at the critical 
point (R c , k c ). The most important fact regarding the EC boundary is that it lies beyond 
the CR boundary for all investigated (L, ^-combinations. Therefore, the EC instability 
occurs only where the rolls are already unstable to CR perturbations. For L = 1, which 
is a typical value for gas mixtures, the CR boundary is still connected to the critical 
point (R c , k c ). However, for smaller L and sufficiently large if) the CR boundary detaches 
from the neutral curve. In this case rolls are no longer a stable form of convection at the 
onset. We observe d square convection at the onset i nstead as it has als o been seen in the 
clear fluid case by Dominguez-Lerma et al. ( 1995 ). Huke et al. ( 2000l ). Miiller fc Liicke 



(|l988t) and lMoses k Steinberd (|i99ll) 



The whole (R, fc)-regions of stable rolls are shown in Figure [7] Inside the closed curve 
defined by the CR boundary rolls are stable to CR perturbations. For large if) and small 
L, the CR boundary shifts to larger wavenumbers and is the only stability boundary 
limiting the region of stable rolls (see L = 0.1 and ip — 0.4). Otherwise, the CR-stable 
region is intersected by the ZZ boundary on the small-fc side. For the presented parameter 
range, we have detected no other instability mechanisms inside the ZZ- and CR-stable 
region. Thus, the regime of stable rolls is limited exclusively by the ZZ and the CR 
boundary in analogy to the pure fluid. 

The CR boundary defines the lowest Rayleigh-Darcy number Ri as well as the highest 
one Rh, for which rolls are stable. As seen from Figure [8j Rh decreases when the sepa- 
ration ratio grows, whereas Ri decreases with if) as soon as the CR boundary becomes 
detached from the neutral curve. Thus, the region of stable rolls shrinks with increasing 
Soret effect. Note that this is also valid where the CR boundary is not detached from 
the ground state and Ri is given by R c . Considering the curvature of Ri(ip) an d Rh(ip) 
it seems that CR-stable rolls cannot exist for too large separation ratios. This conjecture 
is also supported by further numerical tests: We have found no CR-stable rolls for the 
parameters L = 0.1 and if) = 0.5 on an equidistant grid in the (R, fc)-plane with AR = 5 
and Afc = 0.05. 

The wavenumber ki (kh) of the marginally stable rolls at Ri {Rh) is shown in the 
second row of Figure [5] as solid (dotted) lines. Note that ki is equal to k c when the CR 
boundary is connected to the neutral curve. By increasing ip, kh decreases whereas ki 
grows when the CR boundary is detached from the neutral curve. The same trend holds 
for the wavenumber bi (bh) in y-direction of the most critical CR perturbation at (i?;, k{) 
((Rh, kh)), which are displayed as solid (dotted) lines in the third row of Figure [5] While 
bi lies near to ki and bh is similar to kh, they are not exactly the same. 

4.3. Stability boundaries for binary mixtures with small L 
Next, we discuss the choice L = 0.01, which is a typical value for many kinds of liquids. 
Exemplarily for if) = 0.01, the stability balloon is shown in Figure[9] Whereas the EC, CR 
and ZZ boundaries remain qualitatively the same as for larger L, a new stability boundary 
marked by a dash-double-dotted line preceeds the CR boundary for sufficiently large R. 
The corresponding eigenvalue reaches its maximum for d = but for different b in each 
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point of the stability boundary. Moreover, the corresponding perturbations fall into the 
subclass of G-perturbations and are even in x-direction. As we can see, this new boundary 
is rather similar to the CR boundary, and it ends on the latter at smaller values of R. 
However, the difference between the CR boundary and the new boundary is that the 
corresponding eigenvalue of the latter is complex, i.e. the boundary corresponds to an 
oscillatory perturbation. 

In order to clarify the relationship between this new boundary and the CR boundary, 
Figure [TU] shows the real parts of the two eigenvalues with the greatest real parts vs. b 
for different wavenumbers k at parameters L = 0.01, ip — 0.01, R — 290, <f>* = 1 and 
<f? = 0.7. For <fi* — 1, both eigenvalues become positive as a complex pair and thus 
the corresponding stability boundary is oscillatory. But if <fi* = 0.7, the eigenvalues split 
up while they are still negative and only the larger one of the two now real eigenvalues 
becomes positive. Thus, in this case, the corresponding stability boundary is station- 
ary. Both the new boundary and the CR boundary should thus be understood as one 
boundary. We call the oscillatory part oscillatory crossroll (OCR) boundary. 

As a time dependent phenomenon, the OCR boundary can depend on <jf and j a - With 
decreasing cf>* , the OCR boundary moves towards the CR boundary and the i?-values 
increase at which the boundary changes from stationary to oscillatory until the boundary 
is stationary at all R for sufficiently small 4>* . When j a is varied between 10~ 4 and 0, the 
stability boundaries change negligibly. For example, at fixed i?, the wavenumbers change 
less than 0.1%. Moreover, at fixed (R, k), the real parts of the eigenvalues grow less than 
1% with 7„ . The minor influence of j a can be explained as follows: In nondimensional 
units, the diffusion time of the temperature and the concentration field are 1 and 1/L = 
100, whereas the time scale 7 a of the velocity field is at least by a factor 10 4 smaller and 
thus can be neglected. The case that j a has the same order of magnitude as the diffusion 
time of the temperature or the concentration field is rather unusual and lies beyond the 
scope of this work. 

Further research turned out to be difficult. For larger values of ip the boundaries did 
not converge when increasing the number of modes in our model up to N = 40, the 
largest model we used. For large values of (f>* we find the OCR boundary in a large R- 
range while for small values of (jf the CR boundary seems not to be a closed curve any 
longer. 

5. Conclusion 

We studied theoretically roll convection of pure fluids and binary mixtures in the 
Rayleigh-Benard system with porous medium using a Galerkin method. The Soret effect 
was taken into account whereas we restricted the investigations to mixtures of positive 
separation ratio. 

The state of convection has been investigated qualitatively in terms of the streamlines 
as well as the structure of the temperature and concentration field and quantitatively 
in terms of the Nusselt number and the mixing number. Whereas the behaviour of the 
Nusselt and mixing number as well as the behaviour of the temperature and concentration 
field resembles the situation in the clear fluid case, the behaviour of the streamlines 
was found to be different. In the Soret regime, the streamlines were deformed since 
the concentration field transferred its anharmonicity via the buoyancy term into the 
velocity field. For the clear fluid, this deformation is smoothed by the V 2 -operator in the 
diffusion term. However, for a fluid in a porous medium, the diffusion term is replaced 
by a relaxation term in the momentum balance equation deduced from Darcy's law such 
that the deformation persists. At the transition to the Rayleigh regime, the deformation 
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diminishes since form there on the buoyancy force is generated rather by the harmonic 
temperature field than by the anharmonic concentration field. 

Also, we investigated the stability of rolls against arbitrary infinitesimal perturbations. 
The stability region for a pure fluid was found to be restricted by the ZZ instabil i ty an d 
the CR instability only which i s in q ualitative agreement with the one of lStraud (|l974l) . 
However, compared to IStraud (|l974h . our stability balloon is more restrictive on the 
small-fc side and the highest R for which rolls should be stable against infinitesimal 
perturbation is calculated to about 342 instead of 380 stated by Straus. 

In a next step, we extended the stability analysis to binary mixtures with -0 > 0. For 
not too small Lewis numbers {L ^ 0.1), the ZZ and the CR instability mechanism were 
found again to be the only relevant mechanisms. The corresponding stability boundaries 
changed similarly to the clear fluid case: In the Rayleigh regime, the ZZ boundary was 
shifted slightly to smaller k when we increased tp and bent towards smaller k at the 
transition to the Soret regime. The CR boundary was found to be detached from the 
neutral curve for sufficiently large ip. Moreover, the CR stable region shrinked for growing 
separation ratio. Thus, for example, we found no stable rolls for ip 0.5 and L = 0.1. 

For small Lewis numbers such as L — 0.01, we observed a change in the behaviour of 
the CR boundary. For strong thermal driving, the critical eigenvalue corresponding to 
the CR boundary becomes complex. Therefore, the boundary changes to an oscillatory 
one. Up to now, this phenomenon has not been observed in the clear fluid system. The 
point at which the boundary changes from stationary to oscillatory is dependent of ip and 
<p* . When (f>* becomes small enough, the CR boundary is stationary for all values of R. 
Further research for larger tp turned out to be difficult because the computed boundaries 
did not converge. For tp = 0.1 we found that the CR boundary becomes oscillatory in 
a larger i?-range and that the stationary CR boundary does not seem to form a closed 
curve any longer. 

Motivated by the work of Vadasz (for example IVadasz fc Olek (1999)), we retained 
the time derivative in the momentum balance provided with the correction factor j a . 
Whereas, this term was essential in the stability analyses of Vadasz, we did not observe 
any noteworthy influence on our results as j a was varied in the range from to 10~ 4 . 
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Figure 1. In the upper (lower) part, the critical Rayleigh-Darcy number of the ground state 
(corresponding critical wavenumber) is plotted against the separation ratio ip for several L. If the 
thermal driving is weaker than the critical value, the conductive state is stable against arbitrary 
infinitesimal perturbations, otherwise at least one unstable mode exist. 




Figure 2. Structure of roll convection. Shown are from top to bottom the fields of temperature 
and concentration, the streamlines as well as the lateral profiles at midheight, z — 0. Parameters 
are ip = 0.3, L = 0.01, k = n. 
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Figure 3. Nusselt number Nu versus Rayleigh-Darcy number R. The parameters are L — 0.5 
in the top figure and ip = 0.1 in the bottom one. In each case the wavenumber is k = it. PF 
denotes the pure fluid. 
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Figure 5. Upper part: Stability boundaries in the (fc,i?)-plane for the pure fluid. The Eckhaus 
boundary is denoted by the dotted curve, the zigzag boundary by the dash-dotted one, the 
crossroll boundary by the solid one and the neutral curve by the dashed one. The S marks the 
region of stable rolls. Lower part: b ma x on the crossroll boundary. 
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Figure 6. Stability boundaries in the (fc,_R)-plane for several combinations of L and ij) below 
R — 80. The EC boundary is denoted by the dotted curve, the ZZ boundary by the dash-dotted 
one, the CR boundary by the solid one and the neutral curve by the dashed one. The S marks 
the region of stable rolls. 




Figure 7. Full stability balloons in the (k, J?)-plane for several combinations of L and ?/>. The 

curve style is the same as in Figure [6] 
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Figure 8. First row: For three different Lewis numbers, Rh, Ri and R c denoted by the dotted, 
solid, dashed lines, respectively, are plotted versus ip. Second row: The corresponding wavenum- 
bers k denoted in analogy to the first row. Third row: ip versus bh , h and k c denoted in analogy 
to the first row. For further details see text. 




Figure 9. Full stability balloons in the (k, _R)-plane for L — 0.01 and ip = 0.01. The curve style 
is the same as in Figure [6] except that the new OCR boundary is denoted by the dashed-dou- 
ble-dotted line. 




Figure 10. Eigenvalues of the OCR and CR boundary at L = 0.01, ip = 0.01, R = 290, 0* = 1 
or cf>* = 0.7 for different values of the wavenumber k as indicated. Plots show the real parts of 
the two eigenvalues with the greatest real parts. 



